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Abstract The Grad-Shafranov reconstruction is a method of estimating the ori- 
entation (invariant axis) and cross-section of magnetic flux ropes using the data 
from a single spacecraft. It can be applied to various magnetic structures such 
as magnetic clouds (MCs) and flux ropes embedded into the magnetopause and 
in the solar wind. We develop a number of improvements of this technique and 
show some examples of the reconstruction procedure of interplanetary coronal 
mass ejections (ICMEs) observed at 1 AU by the STEREO, WIND and ACE 
spacecraft during the minimum following the solar cycle 23. The analysis is 
conducted not only for ideal localized ICME events but also for non-trivial cases 
of magnetic clouds in fast solar wind. The Grad-Shafranov reconstruction gives 
reasonable results for the sample events, although it possesses certain limitations, 
which need to be taken into account during the interpretation of the model 
results. 

Keywords: Coronal Mass Ejections, Interplanetary; Magnetic fields. Interplan- 
etary; Magnetic fields. Models 



1. Introduction 

Interplanetary coronal mass ejections (ICMEs) are the heliospheric manifesta- 
tions of coronal mass elections (CMEs*) at the Sun. Accroding to Richardson and 
Cane' pOlOp 30% of ICMEs oberved near Earth are MCs. Magnetic clouds as 
defined by |Burlaga et at. (jl98ip are interplanetary structures with dimensions of 



the order of 0.25 AU, which can be identified in in-situ spacecraft observations 
as chunks of solar wind with magnetic field stronger than average, smooth mono- 
tonic rotation of the magnetic field through a large angle, low proton temperature 
and low plasma beta; the most complete up-to-date list of MC signatures can 
be found in IZurbuchen and RichardsonI (j2006p . The present concept of MCs 
assumes a fiux rope to be either tied in the Sun forming a magnetic bottle 
configuration or entirely disconnected from the Sun forming a closed loop. This 
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is proved by observations of bidirectional suprathermal electrons at lAU (Gosling 
et al., 1987). The detection of flare associated electrons within MCs supports the 
magnetic connection to the Sun (e.g. |Farrugia et oT] p993"|) ). 

ICMEs are known to cause the strongest magnetospheric disturbances {e.g. 
Huttuncn et al.. I2002|) . MCs can provide strong southward interplanetary mag- 
netic field. To study the geomagnetic efficiency of an MC it is important to 
know its orientation, shape and size. One can get these critical parameters 
using various flux rope modelling techniques, first attempts of which were made 
by [Burlaga] (|1988p . Numerous flux rope models exist today, such as minimum 
variance analysis (MVA), force-free models such as the Lepping model (Lepping, 



Jones, and Burlaga, 1990), cylinder and torus models by Marubashi and Lepping 



(|20 07!). non-force-free eliptical model by 'Hidalgo, Nieves-Chinchilla, and Cid 



(PD0 2,). kinematically distorted model by |Owens, Merkin, and Riley (2005]), etc. 
All these models fit in-situ observations to an assumed structure of the flux 
rope cross-section. The Grad-Shafranov reconstruction (GSR) technique, on the 
contrary, uses spacecraft observations as initial parameters for the reconstruc- 
tion, thus eliminating the necessity of a priori estimation of the MC boundary. 
GSR was originally developed for reconstruction of flux ropes embedded into 
the magnetopause (Hau and Sonnerup, 1999) and later applied to magnetic 
clouds ( |Hu and Sonnerup, 2002 ). An extented version of GSR useful for multiple 
spacecraft observations was derived by IMostl et al.\ (|2008p . 

In this paper we present improvements to the GSR technique, show exam- 
ples of its usage and discuss main constraints of the method. In Section [2] we 
shortly describe the GSR method and our modifications to it. In Section [3] we 
present examples of events reconstructed with the modified GSR technique and 
in Section [5] we discuss and summarize our results. 



2. GSR and improvements 

The detailed description of the GSR method can be found in |Hau and Sonnerup] 



(|1999p and Hu and Sonnerup (|2002p . Here we just outline the general algorithm 
emphasizing its bottlenecks and possible improvements. 

GSR uses a number of assumptions. Magnetic clouds passing the observing 
spacecraft are assumed to be in magnetohydrostatic equilibrium 

Vp^jxB (1) 

The magnetic field is assumed to have translation symmetry with respect to an 
invariant axis direction, i.e. the fiux rope is assumed to have 2i-dimensional 
structure, where the approximation d/dz = can be used. The whole analysis is 
carried in the deHoffmann- Teller (HT) frame, in which the electric field vanishes 
and thus the magnetic structure can be treated as time-stationary: dB/dt = 0. 

For 2iD magnetic structures with the invariant axis along z equation ([T} can 
be given by the Grad-Shafranov equation 

d^A d^A d / Bl\ 



SOLA: Isavnin_2011a_v2_low.tex; 9 August 2011; 2:19; p. 2 



where A is the magnetic vector potential and the magnetic field vector is B = 
[dA/dy, —dA/dx, Bz{A)]. The plasma pressure, the axial magnetic field com- 
ponent and thus the transverse pressure Pt ^ p + B^/2^o are functions of A 
alone. 

The numerical GS solver is implemented using the Taylor expansions: 

A{x,y±Ay)c^A{x,y)+(^^^ i±Ay) + ^ (^^^ (iAy)^ (3) 



B,{x,y±Ay)~B.,ix,y)+(^^] (±Ay) (4) 



Equation ([2]) expresses an implicit Cauchy problem with a numerical solver 
formed by equations E]). So what we have here is an ill-defined problem 
without boundary conditions, which is subject to growth of singularities after 
certain number of steps of this recurrent algorithm. 



Originally in Hau and Sonnerup (jl999p a finite difference approximation was 
used to estimate d^A/dx^, which is very unstable when dealing with noisy data, 
as most experimental data are. 



In this paper we use smooth noise-robust differentiators ( Holoborodko, 2008 ) 
to improve the stability of the algorithm. Standard finite difference schemes, 
such as central differences, lack high-frequency suppression and may result in 
erroneous results when estimating derivatives in case of noisy data. Noise-robust 
filters, on the contrary, guarantee suppression of high-frequency noise. The ex- 
isting numerical scheme of solving the GS equation ^ suggests multiple differ- 
entiation of initial data, so that instabilities caused by numerical differentiation 
grow like an avalanche. Therefore the use of noise-robust filter can suppress, at 
least to some extent, the growth of such singularities. In Figure [T]the numerical 
second derivative of the function f{x) ~ x^ + NL * RANDN, where NL is 
the noise level and RANDN is a random pick from the normal distribution, is 
shown. From such a simple example it is clear that a noise-robust differentiator 
of the second order gives a much more stable result. 

The determination of the invariant axis is a critical point in the whole recon- 
struction procedure. It is based on the assumption of constant transverse pressure 
and constant magnetic vector potential on common magnetic field lines. The 
search for invariant axis is performed by trial and error. For each test direction 
of the axis magnetic field data are projected to the plane perpendicular to the 
axis. The transverse pressure and magnetic potential are calculated in this plane. 
For the best-fit direction of the invariant axis of the flux rope the Pt{A) curve 
forms two coinciding branches. These two branches represent the motion of the 
spacecraft inward and outward of the flux rope, which correspond to the decrease 
and increase of the distance between the spacecraft and the invariant axis, re- 
spectively. Obviously, the point of the PtiA) curve, which connects two branches, 
correspond to the closest approach of the spacecraft to the invariant axis. Results 
of this search are visualized as a residual map. For each test direction the residue 
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Figure 1. Finite differences method vs noise- robust filter for estimation of the second order 
derivative. 



between inward and outward branches is calculated using the equation 



/\max{Pt) - min(Pt)| 



(5) 



Essentially, the residual map is a contour plot of the residue on top of a 
hemisphere, that represents all possible directions of the invariant axis. The 
initial coordinates used in the residual map are usually defined as follows: y 
is the direction of the maximum variance of the magnetic field, constrained 
by y • Vht = 0; X = — Vht, z = x x y. The direction at 0° longitude and 90° 
latitude is x, the direction at 90° longitude and 90° latitude is y and the direction 
at 0° latitude is z. The search for the invariant axis is performed by stepping 
away from z and calculating residue for each trial direction. We obtain the trial 
direction of the invariant axis in original coordinate system (i.e. GSE, RTN, etc.) 
on each step by rotating z by corresponding latitude and longitude angles. 

An example of the residual map is shown in Figure [5J The black cross shows 
the direction of intcrmidiate variance of magnetic field vector. The direction 
with the minimum residue is considered to be the invariant axis direction and is 
denoted by the thick white dot on the residual map. A problem of this method 
is that the residual maps are occasionally saturated with false possible axis 
directions that correspond to short lengths of branches of the Pt {A) curve and 
thus smaller values of the residue. An example of this is given in Figure [3] (left) . 

To eliminate this issue we combine the residue map with a branch length map. 
The final combined residue is calculated as 



7^ 



(6) 



SOLA: Isavnin_2011a_v2_low.tex; 9 August 2011; 2:19; p. 4 



270 



300 



330 




Figure 2. Residual map for the 1995-10-18 WIND event: original (left) and filtered (right). 
White dot marks the GSR-optimized direction of the invariant axis. Black cross marks the 
direction of intermediate variance of the magnetic field vector. 

where L is the length of the branches (in terms of number of data points) and 
N is the number of observational data points. So 2L/N in ([5]) is a fraction of 
the initial data interval where the coincidence of the branches of Pt{A) takes 
place. Obviously, for the optimal direction of the invariant axis 2L/N shows the 
fraction of the initial data interval occupied by the flux rope. The comparison of 
the original residual map and filtered by equation ([5]) is shown on Figure [5] using 
a WIND event observed on 1995-10-18 as an example. We are using the inverse 
residue 1/TZ in the residue maps. Several false minimum residue directions are 
present in Figure [5] (left) (multiple red zones). After filtering with © we end 
up with a much more clear combined residue map (Figure [2] (right)) with the 
minimum residue direction approximately in the center of the red zone. 

3. Examples of reconstructed events 

In this section we analyze several events using the GSR technique. 

IJian et "all (|2006p presented the classification of MCs based on total per- 
pendicular pressure profile. The total perpendicular pressure is defined as the 
sum of the total magnetic pressure (since the magnetic field does not generate 
pressure force parallel to the mangetic field direction) and the thermal pressure 
perpendicular to the magnetic field direction: 



where j represents different sorts of particles. MCs were divided in three groups 
depending on the shape of the pressure profile: Group 1 with a well-determined 
peak, Group 2 with a plateau observed and Group 3 characterised by decreasing 
pressure. The three groups correspond to a small, medium and large impact 
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Figure 3. Magnetic field and plasma data for the 2008-11-07 STA event. Red vertical lines 
show initial time limits of the MC, green vertical lines show limits of the fiux rope as seen in the 
GS reconstructed magnetic field map. The panels show (from top to bottom): magnetic field 
magnitude (black) and magnetic field components in the RTN coordinates (red: Br, green: 
Bt, blue: Bn), plasma bulk flow speed, proton thermal speed, proton density, plasma thermal 
pressure, total perpendicular pressure of plasma, proton temperature and plasma beta. 



parameters (closest approach of the spacecraft to the axis of the flux rope) 
respectively. We are checking this property throughout our sample events. 

For STEREO events we use the RTN coordinates. In the RTN coordinates, 
R points from the center of the Sun through the spacecraft. T is formed by the 
cross product of the solar rotation axis and R, and lies in the solar equatorial 
plane. N is formed by the cross product of R and T and is the projection of the 
solar rotational axis on the plane of the sky. 



3.1. STEREO-A event on 2008-11-07 



The first event we analyze is a relatively well-defined MC. It was observed on 
November 7th 2008 by the STEREO-A spacecraft (Figure 121). This MC shows 
typical signatures in in situ data: smooth (i.e. low variance) magnetic field ro- 
tation, declining velocity profile (caused by expansion) , low proton temperature 
and low plasma beta. 
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Figure 4. Original (left) and filtered (right) residue maps for the 2008-11-07 STA event. 
White dot marks the GSR-optimized direction of the invariant axis. Black cross marks the 
direction of intermidiate variance of magnetic field vector. 
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Figure 5. Pt{A) curve fitted with polynomial of 3rd order and exponential tail (left) and 
the reconstructed magnetic field map for the 2008-11-07 STA event. The projected coordinates 
shown in the magnetic field map are R (cyan), T (magenta), A'' (yellow). 



The MC originated from a CME event that took place on November 2nd 
2008 at 04:00. The deHoffmann- Teller frame speed estimated for this ICME is 
Vht = [348.0; 3.0; —13.6] km/s in the RTN coordinates with the correlation 
coefficient c = 0.998 ( Khrabrov and Sonnerup, 1998 1. According to the residual 
map (Figure U) the invariant direction of the flux rope is ^htw = 4.3° and 
fRTN = 26.1° in the RTN coordinates. The apex of the MC is the part of 
the flux rope that is furthest from the Sun at any moment of the propagation 
through the interplanetary space. It is natural to expect that the estimated 
invariant axis of the flux rope would tend to be almost perpendicular to the 
radial flow of the solar wind when the spacecraft crosses the flux rope close to 
its apex. And in turn, when the spacecraft crosses the flux rope far from its 
apex, i.e. penetrates through one of its legs, the estimated invariant axis tends 
to be parallel to the radial outflow from the Sun. For this particular event the 
RTN longitude of the invariant axis signifies that the spacecraft intersects the 
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MC far from its apex. The spacecraft crossed this flux rope close to its axis with 
the impact parameter of 0.005 AU. This also agrees with the analysis of total 
prependicular pressure profile (Figure according to which this event falls into 
Group 1 in the lJian et aLI (|2006p classification. The corresponding Pt{A) fitting 
curve for this direction is plotted in Figure [5] (left). Red circles show the in situ 
Pt{A) data when the spacecraft moved inward the flux rope, green circles when 
it moved outward. 

The magnetic field map for this event is shown in Figure[S] (right). It is plotted 
in the MC coordinate system. The picture is the cross-section of the fiux rope 
in the plane perpendicular to the invariant axis. The x-axis is determined by 
the spacecraft trajectory projected to this plane, black arrows show spacecraft 
in situ observations of the magnetic field projected on the same plane. The Sun 
is to the right of the picture. In the upper left corner of the map the projected 
RTN coordinate system is shown. Solid black lines are magnetic equipotential 
lines. The thick white dot, the central point of the flux rope, is the point of the 
maximum magnetic potential. The thick white line is the boundary of the flux 
rope, defined as the absolute minimum of the magnetic potential for which two 
branches of the Pt{A) curve still coincide. The area constrained by this boundary 
may be thought of as an area of reliable GS reconstruction, since for the outer 
part of the magnetic field map the fitting curve for Pt{A) is extrapolated. Using 
this boundary it is possible to get the time limits for the spacecraft passage 
through the fiux rope. These time limits are marked with green vertical lines in 
Figure [3] and arc narrower than those that were estimated by visual analysis. 
This is actually an output parameter of the GSR technique. The first estimate 
of the time limits of the reconstructed event by visual analysis of magnetic 
fleld and plasma data may be rough, since, unlike the flux rope fitting models, 
this method is not that sensitive to the choice of the time boundaries for the 
MC measurements. Note, that for such a well-defined MC the minimum variance 
analysis (MVA) gives the direction of the invariant axis very close to the direction 
estimated with GSR. 

3.2. ACE/WIND event on 2004-11-09: fast ICME 

The MC observed on November 9th 2004 at LI is an example of a fast ICME 
with the average speed higher than 750 km/s (Figure [S]). This MC originated 
from a CME observed on 2004-11-07 at 16:55 UT. Since the propagation speed 
of this ICME well exceeded the average speed of the solar wind, it would be 
expected to produce a leading shock. The analyzed MC was, however, associated 
with two leading shocks, detected by WIND at 09:19 UT and 18:25 UT on 
November 9. Here we reconstruct this MC using data samples obtained by the 
ACE and WIND spacecraft. While multispacecraft reconstruction was already 
performed by IMostl et al.\ (j2008p . we aim here on comparison of these two 
separate reconstruction procedures and their stability. 

Using the residual map analysis (Figure [7]) we estimate the invariant axis 
direction to be Ogse = -27.0°, ipcsE = 42.9° for ACE and Ogse = -15.8°, 
'Pgse = 36.5° for WIND. This MC was crossed by the spacraft rather far from 
its apex, but closer than the MC analyzed in previous section. 
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Figure 8. Reconstructed magnetic field map for ACE (left) and WIND (right) for the 
2004-11-09 event. The projected coordinates shown in the magnetic field maps are Xqse 
(cyan), Ygse (magenta), Zgse (yellow). 



According to lJian et a/.l(|2006p the MCs characterised by the decrease of total 
perpendicular pressure in spacecraft data (as seen in Figure[6|) belong to Group 3 
type of events with the flux rope axis relatively far from the spacecraft trajectory, 
which is proved by our reconstruction results (Figure [51). The impact parameter 
for this magnetic cloud is 0.025 AU for ACE and 0.037 AU for WIND. The dis- 
tance between the spacecraft as seen from magnetic field maps is about 0.008 AU 
which is of the order of the mean distance in the Y Z plane in GSE between ACE 
and WIND, which in turn is equal to 0.002 AU during the event. The step of 
reconstruction along the y-axis used for this event is Aj/ = 0.0001 AU, which 
means that the lack of resolution cannot be the reason of this difference. The 
possible causes for this difference are the inequality between the invariant axis 
directions obtained for two spacecraft and the inaccuracy of the GSR method. 

According to Figure |S] this MC embeds only one flux rope though two leading 
shocks were observed. Note also, that the invariant axis directions estimated by 
MVA and GSR techniques differ by ~ 40 - 50° for this event. 

3.3. STEREO-A event on 2009-07-11: ICME followed by a stream interaction 
region 

The MC observed on July 11th 2009 by the STEREO-A spacecraft is an example 
of an ICME followed by a stream interaction region. As seen in Figure [9] the 
plasma pressure has a well distinguished maximum in the middle of the space- 
craft transit through the MC, which indicates a small impact parameter for this 
event (Group 1 type of events according to IJian et aD (|2006p ). Closer to the 
end of the time interval of a smooth magnetic field rotation the plasma pressure 
grows and well exceeds its value in the central part of the flux rope. This increase 
of plasma pressure is caused by fast solar wind pushing the MC from behind. 
This ICME originated from a CME event on July 7th 2009. The calculated speed 
of the deHoffmann-Teller frame was Vht = [304.8; 4.5; -3.9] km/s in the RTN 
coordinates with the correlation coefficient c = 0.999. In such cases the algorithm 
for invariant axis estimation produces a Pt{A) curve similar to that shown in 
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Figure 9. Magnetic field and plasma data measured on STEREO-A for the 2009-07-11 event. 
Red vertical lines show initial time limits of the MC, green vertical lines show limits of the 
flux rope as seen in the GS reconstructed magnetic field map. The notations are the same as 
in Figure [3] 
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Figure 10. Original (left) and filtered (right) residue maps for the 2009-07-11 STA event. 
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Figure 11. The Pt{A) curve smoothed with running average (left) and fitted with polynomial 
of the 2nd order and exponential tail (right) for the 2009-07-11 STA event. 



Figure 12. The reconstructed magnetic field map for the 2009-07-11 STA event. The projected 
coordinates shown in the magnetic field map are B. (cyan), T (magenta), A'' (yellow). 



Figure [TT] (left). The inward and outward branches of Pt{A) perfectly coincide in 
the central part but differ greatly in transverse pressure in the boundary of the 
flux rope. This happens because the fast solar wind causes transverse pressure 
increase in the rear part of the MC, distorting that part, while the front part 
of the MC remains undistorted. The GS reconstruction may be reliable only 
for the central part of the flux rope where both inward and outward parts of 
the magnetic field lines remain uninfluenced by the fast solar wind. This small 
part of the MC is indicated in Figure [T^ as the boundary of the flux rope. The 
estimated orientation of the flux rope is Ortn ~ —51.7° and (Prtn =33.1° in 
RTN. The flux rope was crossed by the spacecraft at a rather sharp angle far 
from the apex. The impact parameter of the intersection is 0.001 AU. For this 
event the invariant axis estimations produced by GSR and MVA techniques differ 




by - 30°. 
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4. Discussion 



We have presented an overview of the Grad-Shafranov reconstruction technique 
and two improvements to the algorithm. The improvements arc fihering of the 
residual map with the branch length and the estimation of the second order 
derivative of the magnetic potential A used throughout the algorithm with a 
robust noisc-frec filter. We have also conducted a GSR analysis of three sample 
ICMEs that have been chosen to differ greatly in their properties: a clear well- 
defined MC, a fast MC and an MC followed by the fast solar wind. From the 
results of the analysis we see that the GSR gives relevant information about a 
particular MC, although like any data analysis method it should be applied with 
its limitations taken into account. 

The magnetic field maps reconstructed for all three events are almost circular 
in shape. This contradicts the MHD modelling according to which a typical flux 
rope is pancake-shaped ( [Riley and Crooker, 2004] ). This difference was pointed 
out by |Riley et a/.|(|2004p . who conducted a blind test comparison of a simulated 
ICME with various flux rope models. One of the reasons for this apparent flaw in 
GSR may be in its assumptions, i.e. the magnctohydrostatic equilibrium and the 
time stationarity of the MC. The ill-defined mathematical problem solved in GSR 
may also be a reason. Since the Cauchy problem for the elliptic partial differential 
equation ^ is solved numerically without boundary conditions, it is possible to 
obtain almost any result for a given distance from the spacecraft trajectory 



by changing the size of the step along the y direction. In Hau and Sonnerup 



(|1999p this step size was chosen to be Ay/ Ax = 0.1 based on benchmark tests. 
In the benchmarking an analytical solution of the equation V^A = e~^^ was 
used. In our benchmark runs we have obtained better results with Ay/ Ax — 
0.05. Determination of the y-step can only be made empirically for a particular 
reconstruction. 

One possible way of estimating the boundaries of the reconstruction, and 
thus the y-step, is to fix the aspect ratio of the cross-section of the flux rope. For 
instance, in elliptical models the aspect ratio of 4-to-l is often used to mimic the 
pancake shape of the flux rope cross-section {e.g. Mulligan and Russel. [!jOOip . 
The other possible way of estimating the aspect ratio is based on coronagraph 
observations of CMEs. Since the flow of plasma in an MC perpendicular to 
the direction of its propagation is very slow (Owens et al., I2006P the angular 
size of the ICME remains constant. lOwensI (|2008p suggested that coronagraph 
measuremants can be used for estimating the angular size of the CME, i.e. the 
CME width, close to the Sun. Knowing the CME width it is possible to estimate 
the aspect ratio of the flux rope cross-section and use it later as a boundary 
constraint for the GSR. 

Most of the flux rope models are able to give reasonable results only for small 
impact parameters. The GSR possesses this disadvantage too. The reason for this 
is the uncertainty of the definition of the fitting curve for Pt{A) (Figure [T2]) . For 
large impact parameters the high transverse pressure part of the Pt (A) curve is 
undefined for initial observational data and needs to be extrapolated. For various 
extrapolation functions the difference between values of transverse pressure for 
the same values of magnetic potential can well exceed orders of magnitude. 
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Figure 13. A sketch showing uncertainties in Pt{A) fitting for large impact parameters. 



One more reason why the results of GSR may be questioned is that this 
method can show existence of a flux rope in an ICME even when there is none, 
which was pointed out by |Hasegawa et al. (|2007p . 

Despite these shortcomings GSR proves itself as a useful tool for estimating 
the critical parameters of MCs. In this paper we have made GS reconstruction for 
three sample ICMEs under very diverse solar wind conditions. In all these events 
the results of GSR agree with Jian's et al. (|2006p classification of MCs. We were 
able to determine the invariant axis directions for these MCs. The axis directions 
were also estimated using minimum variance analysis (MVA). Note, that for the 
first event the direction of the invariant axis is close to the direction obtained 
in GSR (Figure S]), while for the second and the third events the invariant axis 
estimates given by GSR and MVA differ greatly (Figures [7] (top and middle), 
fTU]) . Thus, for well-defined MCs propagating in slow solar wind with a speed 
lower than the average speed of solar wind GSR and MVA techniques produce 
similar results for the invariant axis direction. While for the non-trivial cases of 
fast MCs and MCs distorted by high speed streams of solar wind GSR is capable 
of giving more relevant results. 
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